arXiv:1504.01628v2 [cs.IT] 13 Oct 2015 


1 


Quickest Eigenvalue-Based Spectrum Sensing using 

Random Matrix Theory 

Martijn Arts, Andreas Bollig and Rudolf Mathar 

Institute for Theoretical Information Technology, RWTH Aachen University, D-52074 Aachen, Germany 

E-mail: {arts, bollig, mathar} @ti.rwth-aachen.de 


Abstract —We investigate the potential of quickest detection 
based on the eigenvalues of the sample covariance matrix for 
spectrum sensing applications. A simple phase shift keying (PSK) 
model with additive white Gaussian noise (AWGN), with 1 
primary user (PU) and K secondary users (SUs) is considered. 
Under both detection hypotheses Ho (noise only) and Hi (signal 
+ noise) the eigenvalues of the sample covariance matrix follow 
Wlshart distributions. For the case of K = 2 SUs, we derive 
an analytical formulation of the probability density function 
(PDF) of the maximum-minimum eigenvalue (MME) detector 
under Hi. Utilizing results from the literature under Ho, we 
investigate two detection schemes. First, we calculate the receiver 
operator characteristic (ROC) for MME block detector based 
on analytical results. Second, we Introduce two eigenvalue- 
based quickest detection algorithms: a cumulative sum (CUSUM) 
algorithm, when the signal-to-nolse ratio (SNR) of the PU signal 
is known and an algorithm using the generalized likelihood ratio, 
in case the SNR is unknown. Bounds on the mean time to false- 
alarm rfa and the mean time to detection ra are given for the 
CUSUM algorithm. Numerical simulations illustrate the potential 
advantages of the quickest detection approach over the block 
detection scheme. 


1. Introduction 

M odern telecommunication systems face ever- 
increasing demands on data rates in order to support 
a growing number of mobile applications, e.g., mobile audio 
and video streaming. Since unlicensed wireless frequency 
spectrum is already very rare, it has been proposed to reuse 
parts of the spectrum which are already licensed. This is 
due to the fact that large parts of the spectrum are not in 
use in certain geographical locations or at certain points 
in time. The term dynamic spectrum access groups efforts 
to make use of such spectrum holes in order to increase 
bandwidth, see [1] for an overview. One subgroup of research 
focuses on opportunistic spectrum access, in which the 
unlicensed users, called secondary users (SUs), determine the 
presence of licensed primary users (PUs) and based on these 
findings decide whether to start unlicensed communication 
in an autonomous fashion. Obviously, a major challenge is 
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the reliable detection of transmission opportunities for the 
SUs, such that ideally no or only a negligible amount of 
interference is caused for PUs. 

Detecting spectrum holes is known as spectrum sensing in 
the literature and there have been a number of publications 
on detection schemes based on various signal features, see [2] 
for a fairly recent review. Among those features, the sample 
covariance matrix and particularly its eigenvalues have been 
used in order to design several detectors, see e.g., [3]-[8]. A 
well known one, which will be relevant in this work, is the 
ratio of the maximum to minimum eigenvalue (MME), i.e. the 
standard condition number (SCN) of the sample covariance 
matrix, see [3], [4]. Since the noise power of the receiver is 
contained in both eigenvalues, it is canceled out in the ratio. 
Thus, designing an alarm threshold for this detector is possible 
without knowledge of the noise power. 

When analyzing communication systems, receiver noise is 
often modeled as additive white Gaussian noise (AWGN). If 
the entries of a matrix A are Gaussian distributed, then the 
matrix AA^ is called a Wishart matrix [9]. On the one hand, 
such matrices occur in the analysis of multiple-input multiple- 
output (MIMO) communication systems as the instantaneous 
MIMO correlation matrix. On the other hand, Wishart matrices 
are also encountered when using the sample covariance matrix 
for spectrum sensing detectors, as will be explained in greater 
detail in Section II. The field of Random Matrix Theory (RMT) 
studies said Wishart matrices and in the recent years there 
has been significant progress on the exact distributions of 
their eigenvalues, see e.g., [10] for a summary. Based on the 
joint ordered eigenvalue distribution, also the distribution of 
the SCN can be found, see [11], [12]. This result, among 
other approximating approaches from RMT, has recently been 
applied to spectrum sensing [13]. In this work, we will utilize 
the SCN distribution in a spectrum sensing problem as well. 

The spectrum sensing detectors discussed above are most 
often employed with a fixed sample size. There, a block of 
samples is taken, the test statistic is calculated and subse¬ 
quently a decision is made by comparing the value of the 
test statistic to a predetermined threshold. One disadvantage 
of said block detection schemes lies in the fact that the number 
of samples taken is independent of the detection difficulty 
of the current situation. Eor example, even if a very strong 
PU signal is present, it will take as many samples until a 
decision is reached as originally planned when defining the 
threshold. This may introduce significant delays, which may 
either result in interference for the PUs or reduce the data rate 
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Fig. 1. Quickest detection problem. At the change time hypothesis 'Hi 
becomes true. The detection algorithm gives an alarm at the detection time 
hence the detection delay or time to detection is — tc- 


by shortening the transmission window of the SUs. 

If one seeks to minimize the detection delay, the hypoth¬ 
esis detection problem can be cast within the framework 
of quickest detection. In contrast to block detection, it is 
assumed that it is known which hypothesis is true and that 
an upcoming change is to be detected with minimum delay. 
Let 7^0 identify the situation that the PU is not transmitting 
and only noise is received. Consequently, let "Hi dehne the 
situation that the PU is transmitting and a signal with additive 
noise is received. Consider without loss of generality that 
there is no PU transmitting (hypothesis "Ho is true). At the 
unknown change time tc a PU starts transmitting, such that 
the corresponding hypothesis Hi becomes true and at the 
detection time the detection algorithm raises an alarm, see 
also Figure 1. The objective in a quickest detection problem 
is to minimize the mean detection delay rj, while having a 
mean time to false alarm Tfa that is greater than a predehned 
value. Quickest detection has been applied to the spectrum 
sensing problem and is also known as quickest spectrum 
sensing in the literature. The hrst publications in this held were 
[14], [15] using Gaussian noise & signal and Gaussian noise 
& sinusoidal signal models, respectively. Quickest spectrum 
sensing using cyclostationary signal features was investigated 
in [16]. Collaboration between multiple SUs using quickest 
detection has been studied in [17], [18]. Also, the case of an 
SU having multiple antennas was analyzed in [19]. 

To the best of our knowledge, so far there has been no 
investigation whether it is feasible to use quickest spectrum 
sensing, while using a test statistic that is based on the 
sample covariance matrix or its eigenvalues. In this work, 
we use a very basic model using two receivers that can 
either be interpreted as one SU with K receiving antennas 
or as K collaborating SUs in close proximity. That is, they 
are assumed to have the same signal-to-noise ratio (SNR). 
It is furthermore assumed, that only one PU is potentially 
present using phase shift keying (PSK) signal modulation. 
With this simple model, we study eigenvalue-based quickest 
spectrum sensing on the basis of the MME test statistic 
mentioned above. The main contributions and the outline of 
this paper can be summarized as follows. First, we formalize 
the signal model and specify the Wishart matrices involved 
under both hypotheses in Section II. Second, we develop the 
exact probability distribution function (PDFs) of the MME test 
statistic under hypothesis Hi for the special case of 2 SUs 
using results from RMT in Section III. Eor the hypothesis Hq 


we recall results known from literature. We also give details 
on the numerical evaluation of the PDEs in the same Section. 
Then, we apply these results to predict the performance of the 
MME block detection algorithm by theoretical calculation of 
the receiver operator characteristic (ROC) in Section IV. In 
Section V we study the quickest detection problem using the 
MME test statistic with and without knowledge of the SNR. 
Eor the hrst case theoretical performance bounds are given. A 
numerical performance evaluation of the MME based quickest 
spectrum sensing detector is subsequently given in Section VI. 
We conclude the paper in Section VII and discuss some future 
research directions. 

II. System Model 

In this Section, we introduce the system model and hx 
notation. We introduce the model in general with K SUs. Erom 
Section III on it will be conhned to two SUs (K = 2) only. 

The basic hypothesis testing problem of detecting the pres¬ 
ence of a PU with K collaborating SUs can be stated in the 
following form: 

Ho : y{t) = w{t) 

Hi : y{t) = x{t) + w{t) , 

where y(f) is a AT x 1 vector of complex baseband samples 
collected by the SUs at time index t G N = {1,2,3,...}. 
Eurthermore, the complex K xl vectors x{t) and w{t) stand 
for the PU’s signal and additive noise, respectively. If Ho is 
true the PU is not transmitting and the SUs receive noise only. 
If Hi is true, the PU is transmitting such that the SUs receive 
a noisy version of the PU signal. Combining N samples into 
a K X N sample matrix yields Y — (y(l),y(2),.. . ,y{N)). 
Similarly, a signal matrix X and a noise matrix W can be 
constructed, so it follows Y = W under Ho and Y — X+W 
under Hi- Using this notation, we can calculate the sample 
covariance matrix as Ry = j^YY". 

In this work we consider a phase shift keying (PSK) signal 
subject to basic additive white Gaussian noise (AWGN). Thus, 
each column of the signal matrix X is assumed to have 
the form Xj = x{j) = Sj Ik for j = 1,... ,N. Here, 
a is the signal-to-noise ratio (SNR), Sj G C is a complex 
PSK symbol on the unit circle (i.e. |sjj = 1) and Ik is a 
column vector of dimension K containing only ones. The 
PSK symbols sj sent by the PU are assumed to be i.i.d. 
and are drawn from a uniform distribution over an arbitrary 
PSK alphabet. The noise matrix is assumed to follow a jointly 
complex circularly symmetric Gaussian distribution, such that 
IV ~ CM{QKy.N ,I K n), where Ik is the identity matrix 
of dimension K. This means, that the real- and imaginary 
part of each matrix entry of W is independently Gaussian 
Af(0,1/2) distributed and the matrix entries are mutually 
independent. Eurthermore, we assume that the PU signal is 
independent of the noise. Under this model, every receiver is 
assumed to have the same SNR. This will be essential for the 
analytical derivation of the PDE under "Hi in Section III. 

The test statistic of the well known MME detector is: 

T=^, (1) 
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where Ai > • • • > Xk are the ordered eigenvalues of the 
sample covariance matrix Ry [3], [4]. Since scaling of the 
sample covariance matrix results in the same scaling of its 
eigenvalues, the ratio is not affected by it. Thus, we will omit 
the normalization factor and use the scaled sample covariance 
matrix R — YY^ in the following. Similarly, the test statistic 
T is only dependent on the SNR a, but independent of the 
actual noise power under both hypotheses. Hence, our system 
model uses the SNR as a parameter directly, as knowledge of 
the actual noise power is unnecessary. 

Under hypothesis Ho the sample covariance matrix is sim¬ 
ply Rq = WW". This is a complex uncorrelated central 
Wishart matrix of dimension K with N degrees of freedom 
[9], which we denote as Rq ~ CWk{N)- 

Under hypothesis "Hi, the sample covariance matrix is 

Ri = {X + W)(X + W)^ 

= XX" + XUU" + UUX" + WVU" . (2) 

The first part of (2) written using the columns Xj of X gives 


III. MME Distributions for Ho & Hi 

As a special case, we have investigated the scenario of two 
cooperating users in more detail, so for the following K = 2. 
For this case, the sample covariance matrices are distributed 
according to Rq ~ CW 2 {N) and i?i ~ C'W 2 {N,aN I 2 ID’ 
respectively. It turns out, that for this simplification the PDFs 
under both hypotheses can be found explicitly. 


A. PDF under hypothesis Ho 

Under hypothesis Ho, the PDF of the test statistic T for 
K = 2 has been found in [11] as: 


MT) 


T{2N) f 

r{N)r{N-i)\ t) \t) 
{N - l)r(2X) (T - 1)2 

[r(X)]2 (T-f i)2tv ’ 
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(4) 


for T > 1. 


N 

N 

= ■ (3) 



Evidently, (3) is constant. Similarly, the middle parts of (2) 
can be expanded as 

N N 

XVF" = H 


WX" = ^ , 

1=1 1=1 

respectively. Since the entries of the vectors Wj are complex 
circularly symmetric Gaussian, their distribution is unaffected 
by complex rotations, which occur by multiplication with the 
unit circle PSK symbols Sj. Therefore the distribution of the 
middle part of (2), XTU" -I- TUX", is identically to the 
distribution of Hence, it makes 

no difference for the distribution of the sample covariance 
matrix (and thereby also no difference for the distribution of 
its eigenvalues) whether an arbitrary PSK signal (including 
4QAM for that matter) or a constant signal is observed. Thus, 
Ri is a complex uncorrelated non-central Wishart matrix of 
dimension K with N degrees of freedom, c.f. [9]. In short 
notation we write Ri ~ CWk{X, n). Here, the parameter 
n = E[(X + TU)]E[(X + TU)]" = is called 

the non-centrality matrix. 

In order to determine the PDF of the test statistic T under 
hypothesis Hi 'm Section III, we need the ordered eigenvalues 
oji > ■ ■ ■ > ujk of the non-centrality matrix TJ. For the rank 
one matrix Jl they are readily given as oji = aKN, while 
UJ2 = ■ ■ ■ = OJk = 0. 


B. PDF under hypothesis Hi 

To develop the PDF of the test statistic T for K = 2 under 
hypothesis Hi, we follow the same way as [20], where a 
version of the PDF was given for X = 2 and generalize this 
result to arbitrary N. We collect the ordered eigenvalues of 
Ri and Ft in the vectors A = (Ai,A 2 )^ and lo = {uji,uj 2 )'^, 
respectively. Then, we can find the joint PDF of the ordered 
eigenvalues of Ri in [12] as: 

2 

A^A) =X„„ |Ui(A)| |F(A,u;)| J]AA,) 

i=i 

= X,„e-(^i+^^)(Ai-A2)(AiA2)('^-') l-F(A,a;)|, 


where |Vi(A)| is the determinant of a Vandermonde matrix 
built from A, |F(A,a;)| is the determinant of a matrix where 
the entry of the i-th row and j-th column follows standard 
generalized hypergeometric functions oiFi{N — l;XjUJi) (see 
[21], (9.14.1)), ^(AA = and 

= [(X-2)!]2(a;i-cc2) ’ 

We are interested in the PDF of the test statistic T, however, 
which is the ratio between the two ordered eigenvalues. By 
substituting Ai = TX 2 we can obtain the desired marginal PDF 
by applying the transformation in (5). Explicit calculation of 
the remaining determinant of (5) yields (6). The hypergeomet¬ 
ric function oHi{a + 1; b) can also be written in terms of the 
a-th order modified Bessel function of the first kind Xa(b), 
c.f. [21], as o^i{a + lA) = o.\ b~^ Ia{2'/b). Substituting (6) 
into (5), using the identity for the hypergeometric function 
and simplifying gives (7). Backing an analytical solution of 
the integral in (7), we use the series expansion of the a-th 
order modified Bessel function of the first kind To( ): 
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{T3 _ 1 ) r(j + 27V - 1) (2a7V)(^-i) 
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( 10 ) 


( 11 ) 


Substituting the Bessel functions by (8) in (7), using that 
iT,T=o^ 3 )iT,Zobi) = and simplifying we 

arrive at (9). The remaining integral in (9) can be found in 
[21] (3.351.3) as 
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a(^+i) 


r((> + i) 

a(6+i) 


to obtain /i(T) for T > 1 and wi ^ UJ 2 in (10). There, we use 
the identity r(a) = (a — 1)! for a G N. As already discussed 
above, the non-centrality matrix O of the PSK signal has rank 
one and therefore wi = 2aN and 0^2 = 0. Setting W 2 = 0 and 
Wi = 2aN in (10), we gain a simpler version of the PDF of 
the test statistic T under hypothesis T-Li for T > 1 in (11), 
which has the SNR a as a parameter directly. 

The PDF /o(T) must be a special case of the PDF f’^{T) 
in the limit for vanishing SNR. Examining (11) in the limit 
for a —?► 0, the only term in the sum which is unequal to zero 
is when j = 1, and it follows: 


lim/“(T) 

a —>-0 


(TV-i)r(2Af) (r-i)2r(^-2) 
[r(Af)]2 (T +1)2^ 


/o(T) . 
( 12 ) 


In Figure 2 the PDFs fo{T) and fi{T) are visualized, the 
latter for different values of the SNR a. 


C. Numerical Evaluation of the PDFs 

Numerical evaluation of the PDFs fo(T) and /f (T) is not 
straightforward, especially for large numbers of samples TV, 
which are relevant in a spectrum sensing context. First, due 
to the factorials and exponents the numerical range needed 
quickly exceeds the range of the IEEE double precision format. 
Second, /“(T) involves an infinite sum and must therefore be 
approximated, e.g., by stopping the summation after a certain 
number of summands. One possibility is to utilize arbitrary- 
precision arithmetic to bypass numerical issues, however, the 
computation time is often unacceptably long. Thus, we will 
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Fig. 2. Plot of foiT) (noise only) and ff(T) for different values of the 
SNR « and At = 500. 


give a brief description of our implementation of the PDEs 
and compare them to empirical PDEs obtained by monte-carlo 
simulations. 

The log-gamma function log(r(a:)) is often used for com¬ 
putations involving factorials or the gamma function itself. We 
used MATLAB for our numerical calculations, where an im¬ 
plementation of the log-gamma function is already provided. 
Therefore, we reformulate fo{T) to perform the numerically 
critical computations inside the exponential function: 

fo{T) = (TV - 1) (T - 1)2 exp [(TV - 2) log(T) 

-2Wlog(T-f l)-flog(r(2Af)) (13) 

-21og(r(TV))] . 

In this way, we can evaluate the function for much higher TV 
than with a direct implementation of (4). 

Since /f (T) contains an infinite sum, we approximate it by 
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Stopping the computation after a finite number of summands. 
Similarly to the approach of fo{T) we reformulate /“(T) to be 
able to use log-gamma functions for numerical computation; 

Js 

/“(T) « (T - 1) ^ exp(/3j -f j log(T)) - exp(/3j), (14) 


/3, = log(rO- + 2N- 1)) - log(r(j + 1)) - {2aN) 
-log(r(7V-l)) + (7V-2)log(T) 

- (j + 2iV - 1) log(T + 1) + (j - 1) log(2a7V) . 

In [20], where a version of /“(T) was developed with 
TV = 2, it was reported that the number of summands can 
be chosen such that a further increase would only add a 
negligible amount (say < 0.5 %) to the summation. However, 
for higher N this criterion fails since the main contribution 
to the summation does not take place in the first summands 
anymore. To achieve a reasonable approximation, we check 
whether the integral over /“(T) is very close to one. In 
Table I, we summarize values for Jg of (14), such that the 
integral deviates less than 10“® from one (using a trapezoidal 
quadrature rule with a bin width of 0.001). 

TABLE I 

Upper bound of summation ( Js) to achieve a reasonable 

APPROXIMATION OF /[* USING (14) FOR DIFFERENT VALUES OF THE 
NUMBER OF SAMPLES N AND THE SNR Ct. 


N l\a 

-20 dB 

-15 dB 

-10 dB 

-5 dB 

0 dB 

50 

10 

20 

30 

60 

200 

too 

20 

30 

50 

200 

300 

500 

30 

60 

200 

400 

2000 

1000 

50 

200 

300 

800 

3000 

5000 

200 

400 

2000 

4000 

20000 

10000 

300 

800 

3000 

7000 

30000 

50000 

2000 

4000 

20000 

40000 

200000 


IV. Application; Block Detection 

In this Section, we apply the results from Section III to 
a block detection scenario, before we turn to the quickest 
detection problem in the next Section. For a block detection 
scheme using the MME detector, the number of samples N 
used in calculating the sample covariance matrix and the 
decision threshold h need to be determined. Intuitively, higher 
N result in less overlap between the PDFs fo{T) and /“(T). 
Thus, choosing N presents a trade-off between the detection 
performance and the time to reach a decision. The threshold 
h, poses a trade-off between the probability of detection Pj 
and the probability of false alarm 

Pf,{h) = 1 - Foih) = 1 - foiT) dT , 

P,{h) = 1 - F^{h) = 1 - fi{T) dT . (15) 

As mentioned above, test statistic T of the MME detector is 
not directly dependent on the actual noise and signal powers, 
but rather on the SNR a. Typically, the SNR a is assumed 
to be unknown beforehand. The decision threshold h is then 
chosen, such that the detector has false alarm rate Pfa. 

Note, that it is assumed in block detection that no change 
of the hypothesis happens within a block. Consequently, when 
the CDEs under hypotheses T-Lq and T-Li are both known, 
the receiver operator characteristic (ROC) can be calculated 
theoretically. Thereby, the design of a detector is simplified 
significantly, as the detection performance can be determined 
beforehand without the need for simulations or empirical tests. 

The CDE under hypothesis T-Lq, i.e. Fq(T) was given in [12] 
as; 

Po(r)=iV„e(p(T)-p(l)), 

where T > 1, 


To verify our implementation including our approximations, 
the PDEs foiT) and /i(r) are plotted in Eigure 3 by eval¬ 
uating (13) and (14), respectively. There, the values of the 
empirical PDEs obtained by digital simulation are drawn in 
circles. 



Fig. 3. Plot of foiT) (noise only) and /f (T) for a = —15 dB and 
N = 500. Circles indicate values taken from an empirical PDF obtained by 
simulation. 


p{T) = Ai{N, N - 1, T) - 2AiiN - 1, N,T) 
+ AiiN-2,N+l,T) 


and 


Ai(a, 6,c) 


( 6 - 1 )! 



(&- 1 ) 

E 


1=0 


ia + j)\c^ 
j\ (c -f l)“+l+i 


We are not aware of an analytical solution to the integral 
in (15). However, with the help of the PDE /“(T) from (10) 
we can numerically evaluate the integral to get Pa- Plotting 
Pd over Pfa for various thresholds results in the ROC. As 
an example, the ROC for N = 500 with various SNRs a is 
shown in Eigure 4. There, results from a digital simulation are 
indicated by circles to verify the theoretical ROC. 


V. Application; Quickest Detection 

This Section applies the results from Section III to introduce 
quickest detection (QD) based on the eigenvalues of the 
sample covariance matrix. In contrast to block detection, where 
the objective is to determine which hypothesis is true, in QD it 
is assumed that it is known which hypothesis is true and that a 
change to the other hypothesis shall be detected with as little 
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Fig. 4. ROC of the MME block detection scheme with N = 500 and various 
SNRs a. Circles indicate values obtained from digital simulation. Compare 
also the corresponding PDFs in Figure 2. 


delay as possible. Contrary to evaluating the probabilities of 
detection Pj and false alarm Pfa as done in block detection, 
the design criterions for QD algorithms are the mean time 
to detection and the mean time to false alarm Tfa. A well 
known algorithm used in QD is the cumulative sum (CUSUM) 
algorithm [22]. It is the optimal QD algorithm, in the sense 
that it achieves minimal rj while satisfying Tfa > c, for any c 
[23]. Without loss of generality, we assume in the following 
that before the change time Q hypothesis Ho is true and at time 
tc transmission of the PU starts, i.e., hypothesis Hi becomes 
true (see also Figure 1). This can be formalized as follows. 

Let Zt be a sequence of random variables, which for each 
t is independently distributed according to the PDF foa{z) 
with parameter Oq before the change time, i.e. for t < and 
according to the PDF fg^ (z) with parameter 0i at and after the 
change time, i.e. for t > Q- Samples z{t) are taken in order 
to detect the change from Ho Hi. The CUSUM algorithm 
relies on the log-likelihood ratio of the received sample at 
time index t: l^{t) = log {fe^{z{t))/ fgg{z{t))). Over time, 
lz{t) shows a positive drift if Hi is true and a negative drift 
if Ho is true. A cumulative sum of the log-likelihood ratios 
of consecutive samples is formed by the CUSUM algorithm. 
However, only those consecutive samples are used for the 
summation, which result in the largest sum. In other words, 
only the time interval showing a consistent positive drift is 
considered for summation. A convenient recursive formulation 
of the CUSUM algorithm can be given as [22], [24]: 


9z{t) 


t 

max y log I 

0<m<t ^ ' 

[gzit -1) + izit)] 


feJADY] 

feo{z{j))J 


(16) 


where the initial value is defined as Pz(0) = 0 and the positive 
part is defined as = max(-,0). At each time index t 
the value gz{t) is compared to a threshold h > 0 and if 
gz{t) > h the algorithm decides that a change to Hi has 
happened. As can be seen from (16), gz{t) may never reach 
negative values. Since lz{t) shows a negative drift under Ho 
and we are interested in detecting a potential change from 
Ho —> Hi, accumulating negative values would result in a 


much longer time to reach the threshold after the change time 
tc and thus in a larger mean time to detection rj. Note, that 
the parameters 9o and 9i must be known to use the CUSUM 
algorithm. If the parameters are unknown, the algorithm can 
be adapted to use the generalized log-likelihood ratio. We will 
study the application of QD using the MME test statistic with 
and without knowledge of the SNR in the following. 


A. Known SNR 


First, we study the model from Section II for the case that 
the SNR a is known to the receiver. To adapt the results from 
Section III, we define a time-dependent version of the test 
statistic T from (1). Introducing a block index fc G N, we 
define A(fc) as the vector of ordered eigenvalues of the scaled 
sample covariance matrix of block k: R{k) = Y{k)Y^(k), 
where Y (k) = (y((A: — 1)W 4-1), ••• ,y{kN)). In other 
words, A(fc) contains the eigenvalues of the fc-th consecutive 
sample covariance matrix built from a block of N non¬ 
overlapping samples. Thus, since we assume K = 2, the test 
statistic follows as: 


T{k) = 


Ai(fc) 

A2(fc) ■ 


To take over existing results from QD, we make an addi¬ 
tional assumption, which is similarly present in block detec¬ 
tion: Q € {(A:—1)A^-|-1 I k G N}. This means that no hypothe¬ 
sis change may happen within a block. Depending on whether 
Ho or Hi is true, T{k) is distributed according to fo{T) or 
fi{T), respectively. As shown in (12), lim /“(T) = fo{T) 

CK—>-0 

so the only parameter influenced by the change is the SNR a. 

The log-likelihood ratio can be found by inserting (4) and 
(10) into its definition and simplifying as seen in (17). There, 
we use the Pochhammer symbol defined as (a)j^ = {a + b — 
l)!/(a — 1)!. Note, that (17) can be evaluated numerically 
using a similar approach to the one used for the PDFs in 
Section III-C. Using the log-likelihood ratio, we can give the 
CUSUM algorithm for our model as: 


gik) = [gik-l)+lik)]+ . (18) 


When designing the CUSUM detector, a method for finding 
a suitable threshold h is desirable. Since the exact calculation 
of the mean time to detection Td and the mean time to false 
alarm Tfa is very complicated in general, we turn to bounds 
already studied in the literature. Of special interest are an 
upper bound on Td and a lower bound on Tfa to ease the 
design of the threshold h. Note, that we are evaluating the 
bounds using the timescale of the block index k. To evaluate 
the bounds on a sample based timescale (time index t in 
Section II), the bounds must be multiplied by the number of 
samples N used for calculating the sample covariance matrix 
in each block. 

Starting with the upper bound on Td, we And from [24], 
[25], that 

(h+iim) 


Td < 


E/r im] 


( 19 ) 


where 


lif) = sup Fif [l{k) — S I l{k) > (5 > 0] . 

i5>0 
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l{k) = log 


{ fnTm \ 

\ fo{T{k)) ) 


( “ {nky - 1) (2a7V)0-i) (2iV)(^_,) 

yiTik) - 1) ^ j! (r(fc) + 1)0-1) 


(17) 


Here, we denote the expectation over the PDF /(T) as E/[-]. 
Likewise, a lower bound on Tfa can be found in [24]: 


1 

- Ef, m] 



+ h + 7(/o) 


( 20 ) 


where (/? < 0 is the single non-zero root of E/(,[e 
which can be found by solving: 


^/o(r)dT = l. (21) 

Obviously, (p = —1 solves (21), since only the integral over 
fi{T) remains. Thus, (20) can easily be simplified to: 

Neither [l{k)], E/o [l{k)], 7 (/o) nor 7 (/f) lend them¬ 
selves to proper analytical treatment, but can be calculated 
numerically. 

A second lower bound on Tfa reported in the literature (see 
e.g., [24]) is simply 

7-fa > . (23) 


Although this bound is very simple to compute, we found that 
for our model and the threshold values we considered in the 
numerical evaluation, this bound turns out to be very loose. 
Since (22) is much tighter to the value of Tfa obtained by 
simulation, we will omit (23) in our numerical evaluation. 


B. Unknown SNR 


A more realistic scenario is that the SNR is unknown to 
the receiver beforehand. In this case, the CUSUM algorithm 
cannot be used directly. However, a similar algorithm based 
on the generalized likelihood ratio test (GLRT) can be utilized 
instead. The GLRT allows to perform likelihood ratio tests 
when the parameters of the involved PDFs are unknown, 
by estimating these parameters using a maximum likelihood 
estimation (MLE) beforehand [26]. In our scenario, the only 
unknown parameter is the SNR and the GLRT for this case 
can be given as: 

where a is the estimated SNR. Utilizing the GLRT test on 
an i.i.d. sequence of samples, the generalized likelihood ratio 
(GLR) algorithm can be constructed as [24]: 


9G{k) 


max sup log 

0<m<k ^ 


( A fHnj)) 


k 

= max sup > log 

0<.m<.k ^ 

- - OL 


f fHnm 


(25) 


We are not aware of an analytical form of the supremum 
in (24), therefore we evaluate it numerically. Note also, that 
the GLR algorithm of (25) cannot be written in a recursive 
form. In contrast to the CUSUM algorithm, the GLR algorithm 
requires memory of all previous samples. This is due to the 
numerical evaluation of the supremum and the subsequent 
maximization in (25), which also has a higher computational 
complexity. 

In Figure 5 the mean time to detection obtained from 
digital simulation is plotted over the threshold for the CUSUM 
and the GLR algorithm. Also, the theoretical upper bound on 
T(j from (19) is depicted. Note, that the upper bound is only 
valid for the CUSUM algorithm. Similarly, in Figure 6 the 
mean time to false alarm Tfa obtained by digital simulation 
is shown for both the CUSUM and the GLR algorithm with 
a logarithmic scaling of the ordinate. Again, the theoretical 
lower bound on Tfa from (22), which is only valid for the 
CUSUM algorithm, is evaluated there. The values of and 
Tfa shown in the Figures 5 and 6 are evaluated on a block-wise 
timescale, as already explained above. 



Fig. 5. Theoretical upper bound from (19) on the mean time to detection 
Td of the CUSUM algorithm. Simulated performance of both the CUSUM 
and GLR algorithm is shown for N = 500 and SNR a = —15 dB. 
The GLR numerically evaluates the supremum over d within the interval 
[-20 dB, -5 dB] in 0.1 dB steps. 

As can be seen in both Figures 5 and 6, the GLR algorithm 
shows a more aggressive behavior compared to the CUSUM 
algorithm. This can be explained by two effects, predominantly 
present at the beginning of the algorithm run-time. Firstly, 
the SNR estimate a used by the GLR algorithm (25) may 
over- or underestimate the SNR. If, for instance, the hrst 
couple of sample covariance matrices result in a greater than 
average test statistic T{k) for the respective SNR, the GLR 
will overestimate the SNR a and hence the log-likelihood 
ratio will be greater than it would be for the true SNR 
a. Secondly, when the GLR algorithm (25) determines its 
parameter m, it may also choose larger values for m than 
it would if the SNR estimation d had returned the true SNR 





















Fig. 6. Theoretical lower bound from (22) on the mean time to false alarm 
Tfj, of the CUSUM algorithm. Simulated performance of both the CUSUM 
and GLR algorithm is shown for N = 500 and SNR a = —15 dB. 
The GLR numerically evaluates the supremum over a within the interval 
[—20 dB, —5 dB] in 0.1 dB steps. Note the logarithmic scaling of the ordinate. 

a. In Figure 7, a typical run of both the CUSUM and the 
GLR algorithm are depicted. The bottom of Figure 7 shows 
the internal parameters of the GLR algorithm. Both effects 
can be observed as discussed above. With advancing run¬ 
time, the SNR estimation d converges to the true SNR a. 
Consequently, also the behavior of the GLR algorithm more 
closely resembles the behavior of the CUSUM algorithm. If 
the SNR estimation of the GLR algorithm returns the true 
SNR, both algorithms behave identically, as can be seen by 
comparing (18) and (25). 



block index k 


Fig. 7. Comparison of a typical run of the CUSUM and the GLR algorithm 
with N = 500 and SNR a = —15 dB. The internal parameters a. and m of 
the GLR algorithm are shown in the bottom plot. 

VI. Numerical Evaluation 
A fair comparison between the GLR algorithm, which is 
a quickest detection algorithm and the MME block detection 


scheme is difficult to obtain due to the assumptions present in 
the two different detection problems. In the following evalu¬ 
ation we design the MME block detector with a block length 
of 10® samples and probability of false alarm Pfa = 0.0147. 
Using the results from Section IV, we find the probability 
of detection = 0.9269 for the SNR a = —20 dB. 
When using the GLR algorithm, a run-time has to be decided 
because the algorithm will eventually generate a false alarm. 
Evidently, this run-time is to be chosen such that the mean time 
to false alarm Tfa is significantly higher than the algorithm 
run-time. We performed a monte-carlo simulation with the 
GLR algorithm for 1000 random seeds, with N = 10000 
samples in each block and a variety of threshold values. 
Eurthermore, two versions of the algorithm per monte-carlo 
instance are computed. The first one receives samples under 
hypothesis T-Lq, i.e. it only receives noise. The second one 
receives samples under hypothesis Tfi, so it receives signal 
with additive noise. Eor the detection performance this is a 
worst-case analysis, since the GLR algorithm is initialized with 
the value zero. If the hypothesis would switch from Tfg to "Hi 
during the simulation, it might detect a longer consecutive 
positive drift due to the effects of the noise, resulting in 
potentially quicker detection. To calculate probabilities of false 
alarm and detection for the GLR algorithm, respectively, we 
consider different run-times and calculate how many of the 
1000 monte-carlo instances have raised a false-alarm / have 
correctly detected the signal for a given run-time, threshold h 
and SNR a. Thus, we can compare the detection performance 
of the GLR algorithm to the fixed performance of the MME 
block detector. Eor the numerical evaluation of the supremum 
in (25), an SNR range of [—25, —5] dB was evaluated in 0.1 dB 
steps. In Eigure 8, the GLR performance is depicted for several 
SNRs. The threshold h = 4.5 was chosen, such that at a run¬ 
time of 10® samples, i.e., the block size of the MME detector, 
the probabilities of false alarm are approximately equal (GLR: 
Pfa = 0.0140 vs. MME: Pfa = 0.0147). As expected, the GLR 



Fig. 8. Performance of the GLR algorithm, evaluated as Pd and Pfa over 
the algorithm run-time for the threshold h = 4.5. The gray circular markers 
and the thin solid lines indicate the performance of the MME block detector 
designed for a block length of 10^ samples and Pfa = 0.0147, which results 
in Pd = 0.9269 at SNR o = -20 dB. 
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quickest detection algorithm is able to detect higher SNRs 
much earlier than the MME detector at a lower or comparable 
false-alarm rate. However, at low SNRs the MME detector 
is faster. When considering realistic detection scenarios, the 
PUs SNR cannot be known beforehand. So the MME detector 
is designed to have a specihc performance for the lowest 
SNR that has to be reliably detected. Evidently, in practice 
it would be more precise to consider an SNR interval, say 
e.g., [—5, —20] dB. If no further knowledge is available, we 
may model the SNR with a uniform distribution over said 
interval. Here, we see the advantage of the quickest detection 
approach. Eor a wide range of SNRs, it offers faster detection 
at comparable false-alarm performance, thus wasting less time 
on the detection process. On the one hand, this leads to 
higher efficiency in using transmission opportunities, thereby 
increasing transmission rates. On the other hand, interference 
for the PU may be decreased by diminishing the reaction time 
to free the channel when the PU initiates a communication. As 
can be seen in Eigure 8, there is a gap of approximately 2 dB 
in the low SNR performance of the GLR quickest detection 
algorithm compared to the MME block detector. This is due 
to the fact that the GLR algorithm does not use the entirety of 
samples, but rather blocks of N = 10000 samples to calculate 
the sample covariance matrices. Thus, in choosing the number 
of samples N, one has to trade-off the detection delay for 
higher SNRs and the performance gap compared to the MME 
block detector for low SNRs. However, it must be stressed 
that this comparison is not completely fair, since the GLR 
algorithm introduced in (25) essentially throws data away. This 
stems from the CUSUM analysis, where the bounds rely on 
the assumption that the samples of the test statistic used in 
the CUSUM algorithm are i.i.d.. Eurther improvements in this 
direction are discussed in Section Vll. 

As for all detection algorithms the choice of the detection 
threshold has a huge impact on the detection performance. 
This effect can already be seen in Figures 5 and 6. To visualize 
the performance trade-off between Pj and Pfa directly, three 
different choices of the threshold h are compared in Figure 9. 


VH. Conclusion & Future Work 

In this work, we have investigated the potential for using 
quickest detection based on the sample covariance matrix 
for spectrum sensing applications. First, we have specihed a 
simple phase shift keying (PSK) model with additive white 
Gaussian noise (AWGN), consisting of one primary user 
(PU) and K secondary users (SUs). Second, we characterized 
the Wishart distributions of the eigenvalues of the sample 
covariance matrix under both detection hypotheses Ho (noise 
only) and Hi (PU signal + noise). Third, we derived analytical 
formulations of the PDF of the maximum-minimum eigen¬ 
value (MME) test statistic, which is also called the standard 
condition number (SCN) of a matrix, under both hypotheses 
Ho and Hi for the special case K — 2. These PDFs have the 
number of samples N in the sample covariance matrix and 
the signal-to-noise ratio a as parameters. Then, we discussed 
the numerical computation of said PDFs. We put these results 



GLR run-time .^qS 


Fig. 9. Performance of the GLR algorithm, evaluated as and Pfa over 
the algorithm run-time for different thresholds h at SNR a = —17 dB. 


to application by explicitly calculating the ROC for the well 
known MME block detector without the need of expensive 
simulations. Subsequently, the main application of these re¬ 
sults is discussed: quickest detection based on the sample 
covariance matrix, i.e. eigenvalue-based quickest detection for 
spectrum sensing. Two quickest detection algorithms were 
introduced, which rely on the MME test statistic. The CUSUM 
algorithm is applicable when the SNR of the PU signal is 
known and a GLR algorithm was deduced to cope with the 
situation when the SNR is unknown. Bounds on the mean time 
to false-alarm Tfa and the mean time to detection have been 
provided for the CUSUM algorithm. Numerical simulations 
illustrate the potential advantages of the quickest detection 
approach compared to the block detection scheme for spectrum 
sensing applications. It turns out, that for a wide range of SNRs 
the introduced GLR quickest detection algorithm offers faster 
detection at better or comparable false-alarm performance than 
the MME block detector. 

Future research may focus on hnding similar bounds on 
Td and Tfa for the GLR algorithm. Also, at very low SNRs 
the block detector is still faster than the GLR algorithm. This 
results from the fact that the GLR algorithm calculates its 
sample covariance matrices block-wise and hence uses less 
data than the block detector. One possibility of improvement 
would be to get rid of the block structure introduced in 
Section V and perform updates of the samples covariance 
matrix on a sample basis. However, this would introduce 
stochastic dependencies in the eigenvalues, so that the behavior 
of both algorithms would be changed fundamentally. Thereby 
the presented bounds would be rendered invalid. Since the 
presented GLR algorithm is computationally expensive, it 
would be interesting to investigate efficient heuristics which 
can be used in practical scenarios. Additionally, it would 
be interesting to extend the simple channel model to more 
realistic fading channel models. Finally, it may be relevant to 
examine other test statistics based on the (eigenvalues of) the 
sample covariance matrix for quickest spectrum sensing. 
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